First-principles calculation of the thermodynamics of In^Gai-^N alloys: 

Lattice vibrations effects 
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The thermodynamics properties of the wurtzite and zinc-blende In^Gai-^N alloys are calculated 
using first-principles density-functional calculations. Special quasi-random structures are used to 
describe the disordered alloys, for x = 1/4, 1/2, and 3/4. The effect of lattice vibrations on the phase 
diagram, commonly omitted from semiconductor alloy phase diagram calculations, are included 
through first-principles calculations of phonon spectra. Inclusion of lattice vibrations leads to a large 
reduction in the order- disorder critical temperature (~ 29% and ~ 26% for the wurtzite and zinc- 
blende structures, respectively) and changes the shape of the solubility and spinodal curve through 
changes in the entropies of the competing phases. Neglect of such effect produces significant errors 
in the phase diagrams of complex ordered semiconductor compounds. The critical temperature for 
phase separation is 1654 K (1771 K) for the wurtzite (zinc-blende) structures. The predicted phase 
diagrams are in agreement with experimental measurements on MOCVD In^Gai-^N films. 
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I. INTRODUCTION 



The group-Ill nitrides GaN and InN and their al- 
loys have been receiving considerable attention for high- 
power, high-frequency, and high-temperature optoelec- 
tronic applications, such as light emitting and laser 
diodes. By changing composition one can, in principle, 
continuously tune the band gap from 0.8 eV (InN)i to 
3.5 eV (GaN) — a range that spans much of the visi- 
ble spectrum. However, because of the large lattice mis- 
match between GaN and InN (~ 11% for the a and 
c-directions in the wurtzite structure), the solid solu- 
tion has a tendency to undergo phase separationi2iii4 
GaN, InN and their alloys (In^Gai-^N), grown un- 
der ambient conditions, typically assume the hexagonal 
wurtzite (WZ) crystal structurei^^ It is possible, how- 
ever, to grow thin epitaxial films of In^Gai-^N with 
the cubic zinc-blende (ZB) structure. 7 ! 8 ! 9 ! 10 ! 11 ! 12 ! 13 ! 14 
The phase stability of In^Gai-^N alloys has been 
the subject of several experimental 3 ^ and theoretical 
studieo 15 ! 16 ! 17 ! 18 ! 19 ! 20 ! 21 ! 22 ! 23 ! 24 . For example, there have 
been more than ten calculations of the In^Gai-ajN 
pseudo-binary phase diagram. These have included 
valence-force field (VFF)^*^^^, tight-binding, and den- 
sity functional calculations^, sometimes combined with 
atomic-scale Monte Carlo simulations^. Surprisingly, 
there is considerable discrepancy in the predicted mis- 
cibility gap. The tendency for phase separation is 
largely driven by the relatively large misfit strains^. 
The predicted critical temperature varies over a 2000 K 
range, with most calculations predicting 1500 K < T c < 
2000 K. The source of this discrepancy may be associ- 
ated with the accuracy of the description of the atomic 
interactions, the need for a statistically meaningful de- 
scription of the distribution of Ga and In on the solid so- 



lution cation sites and the contribution to the free energy 
associated with configurational entropy^. Additionally, 
the vibrational contribution to the free energy has never 
been included. In the present paper, we report the re- 
sults of the most rigorous calculation of this phase dia- 
gram performed to date. The results clearly demonstrate 
the importance of the vibrational contribution in deter- 
mining phase diagrams in semiconductor alloys, such as 
In x Gai_ x N. 

In this work, we employ accurate density-functional 
theory (DFT) to examine the thermodynamic properties 
of both wurtzite and zinc-blende In^Gai-xN. We directly 
compare the formation enthalpy and free energy differ- 
ences determined using exactly the same pseudopotcn- 
tials and energy cutoff for both structures. Using these 
results, we predict the binodal and spinodal on the binary 
phase diagram; delineating the miscibility gap and the re- 
gion where the solid solution is unstable. The simulation 
method uses the special quasirandom structure (SQS) 
formalism 2 ^ 2 ^ to faithfully (and efficiently) represent the 
structure of the random alloy. The effect of lattice vibra- 
tions, often neglected in semiconductor alloy thermody- 
namic calculations, is included in a harmonic approxi- 
mation where the dynamical matrix is determined using 
density- functional theory methods. Even though lattice 
vibration effects have been shown to be significant for 
the thermodynamic properties of some intermetallicsS 7 -, 
the effects of vibration have been commonly neglected 
for semiconductor alloys. Such effects were not previ- 
ously included, in part, because an accurate treatment 
of lattice vibrations requires substantial computing re- 
sources, above that needed for modeling the effects of 
alloying on the heat of formation. Advances in first- 
principles computational methods and hardware make 
routine inclusion of vibrational effects in alloy phase di- 
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agram calculations possible. As wc demonstrate below, 
lattice vibrations can drastically alter calculated thermo- 
dynamic properties and phase diagrams of In^Gai-^N. 
This work provides the first direct evaluation of the sig- 
nificance of lattice vibration effects on the phase diagram 
of compound semiconductor alloys — a large and tech- 
nologically important class of complex alloys. Earlier 
work on solid solution semiconductor alloys Si-Ge^ and 
Gai-a-In^P— suggested that lattice vibration effects are 
negligible. On the other hand, inclusion of these effects 
on metallic compounds (intermetallics) have been shown 
to be important in determining phase stability in the Cd- 
Mgffi, Al-So2i, Cu-AuS, and Al-Cu^ systems. However, 
vibration effects have been shown to be relatively unim- 
portant in other metallic systems that exhibit compound 
formation, e.g., Ni-Al2i and Pd-Vi&. The previous work 
on the phase stability of metals suggests that vibrational 
effects may or may not be important in those systems. 
Further, there is no evidence to support the idea that one 
can extrapolate from metals to compound semiconductor 
alloys. 

The solid state In^Gai-^N phase diagrams, deter- 
mined herein, represent the most rigorous determinations 
reported to-date. The resultant phase diagram shows 
that the wurtzite crystal structure is more stable than the 
zinc-blende structure over the entire In 2 .Gai„ a; N compo- 
sition and temperature ranges examined. 



II. PREVIOUS THEORETICAL 
CALCULATIONS 

Due to its low computation cost and relative accuracy, 
valcnce-force-field (VFF) models^ have been used ex- 
tensively to studjii^iSiiLiS the In^Gai-^-N alloy system. 
One of the first In^Gai-^N phase diagram calculations 
was performed by Ho and Stringfellowi^, who employed a 
VFF model to investigate the solid phase miscibility gap 
in zinc-blende In a; Gai_. r N. They found a critical temper- 
ature of 1523 K and that at the growth temperature of 
800° C, the InN solubility in GaN is less than 6%. Using 
interaction constants deduced from first-principles cal- 
culations, Saito and Arakawai& simulated the wurtzite 
In x Gai_ x N alloy system using the VFF model and ob- 
tained a critical temperature of 1690 K. Takayama et 
aiJi, using a modified VFF for In^Gai-^N, found a crit- 
ical temperature of 1668 K for the zinc-blende structure 
and 1967 K for the wurtzite structure. 

Several first-principles calculations have also been per- 
formed. Van Schilfgaardc et alM> relaxed a 32-atom 
structure using a tight-binding method and determined 
the heat of formation using density-functional theory. 
Combining this result with an ideal mixing estimate of 
the configurational entropy led to critical temperatures 
of 3980 K and 2290 K for the wurtzite and zinc-blende 
structures, respectively. These values are considerably 
higher than those obtained from the VFF models. Teles 
et alml investigated strained and relaxed In^Gai^^N al- 



loys using a pseudopotcntial plane-wave approach with 
a generalized quasichemical approximation (GQCA) for 
mixing 37 . They obtained a critical temperature of 1295 K 
for the zinc-blende In^Gai-^N, which is lower than the 
VFF resultsi^ii. Another first-principles calculation 23 
also reported a similar value for the critical temperature 
(i.e., 1400 K) for zinc-blende In^Ga^^N. 

Grosse and Neugebaues 3 ^ assessed the limits and accu- 
racy of the VFF models for In^Gai^N alloys by studying 
various ordered structures using first-principles calcula- 
tions. They found that while VFF works reasonably well 
for some zinc-blende III/V semiconductors, significant er- 
rors occur for group III nitrides (GaN, InN). In partic- 
ular, the formation energies for alloys are significantly 
underestimated (by 14 meV/cation). These findings sug- 
gest that the critical temperature of In^Gai-^N obtained 
from first-principles calculations should be higher than 
those from VFF calculations. However, this conclusion 
is inconsistent with the predictions of VFF^*iS*ii and 
first-principles results^SiS 3 . (i.e., the T c values from VFF 
studies exceed those from first-principles calculations). 
Our results, presented below, shed light on this discrep- 
ancy. 



III. SPECIAL QUASIRANDOM STRUCTURES 

The effects of disorder on the stability and structure 
of alloys have received considerable attention over the 
past few decades. The coherent-potential approximation 
(CPA) 3 ^ for describing such effects is a single-site the- 
ory that treats the random A^Bi-^ alloys by considering 
the average occupations of lattice sites by atoms A and 
B. Hence the effects of local environment such as charge 
transfer, local chemical environment, and local structural 
relaxation are ignored. However, local environment ef- 
fects have been shown to be important in determining the 
thermodynamics and electronic properties of alloys 2 ^ 2 ^. 
One approach to incorporate the local environment effect 
is through first-principles density-functional theory. One 
can, in principle, model a disordered A x J5i^ x system us- 
ing a plane- wave density-functional methods with a large 
supcrccll in which A and B are distributed at random on 
lattice sites. Unfortunately, since computational costs in 
such methods usually scale with the number of atoms in 
the supercell, N as N 3 , one is usually forced to use a 
small supercell. Small supercells invariably lead to poor 
statistics and spurious correlations (i.e., this approach 
leads to a chemical environment about each site which 
is different from that in a large, truly random system). 
Zunger et alm^^ recognized that most physical prop- 
erties are governed by atomic interactions between near 
neighbors and developed a method to generate atomic ar- 
rangements in small supercells that correctly reproduce 
the short and intermediate range correlations that exist 
in random alloys. These structures are called special- 
quasirandom-structurcs (SQS). 

A useful approach to defining correlation functions in 
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crystals with different site occupancies is to assign each 
site i a variable according to the type of occupying atom, 
say Si — +1 for an A atom and s,- = — 1 for an B atom. 
Next, we make a list of pairs of atoms in our structure 
that are separated by k or less nearest neighbor distances 
and evaluate the product tt = SiSj for that particular pair 
of atoms on sites i and j. Averaging this value over all 
such pairs in the structure yields the correlation function 
hb.fc, where the two indicates we are considering pairs. 
We can use a similar procedure for triplets, quadruplets,., 
m-lets of atoms to produce a series of correlation func- 
tions H m ,k- If the alloy were of infinite extent and the 
site occupancies were truly random, all Tl m ,k would be 

U m , k {x) =n m {x) = [(+l)z + (-l)(l-z)] m = (2a-- 1)™, 
where x is the concentration of A atoms in the system. 

We have implemented the SQS algorithm for different 
unit cell sizes and shapes, generating all possible site oc- 
cupancies for the chosen x, calculating II m ,fc and keeping 
the configuration for which the set of n mi jt(x) is closest 
to H m (x). In practice, we use only m = 2 and k < 10. 
We repeat this procedure for different unit cell shapes, 
as described by the supercell lattice vectors ai, &2, and 
a3, and retain the one that leads to the best correlation 
with random occupancy (using our implementation of the 
procedure described in Refs. lioEH) . 

In this work we have focused on wurtzite and zinc- 
blende crystal structures with 16 cations and 16 anions 
per unit cell. The occupancies, in our case correspond 
to Ga or In on the cation sites. The descriptions of 
the hexagonal-close-packed (hep) and face-centered-cubic 
(fee) sublattices are shown in Tables [I] and EI respec- 
tively. For x = 1/4 (or x = 3/4), the wurtzite structure 
delivers exact matches up to 112,5 and up to 112,7 f° r 
x = 1/2. For the zinc-blende structure, we find struc- 
tures that yield correlation functions that exactly match 
the random case up to 112,3 for x = 1/4 (or x — 3/4) 
and to 112,7 for x = 1/2. Unlike in previous studies 
of Ina;Gai_ x N alloys^ 3 -, the present analysis provides a 
quantitative measure of the degree to which the simula- 
tion cell resembles a random alloy. Based on the quoted 
pair statistics, we conclude that the unit cells and occu- 
pancies employed herein provide good representations of 
the set of wurtzite and zinc-blende random alloys. 



IV. FORMATION ENTHALPY CALCULATIONS 

We first optimize the zinc-blende (wurtzite) GaN 
and InN structures with respect to the atomic coor- 
dinates and the unit cell parameter (s) a (and c) us- 
ing the first-principles density-functional 4 ^ 3 , method 
within the local density approximation (LDA). This was 
done using the plane-wave pseudopotential Vienna ab 
initio simulation package (VASP)^non- norm-conserving 
pseudopotentials^ The importan t 38 ! 46 3d states for Ga 
and 4d states for In are included as valence electrons in 
the pseudopotentials. A relatively high cutoff energy of 
31.96 Ry is used throughout this work. A Monkhost- 



Pack mesh of 7 x 7 x 7 is used for the 4-atom unit cell 
wurtzite InN and GaN systems. A 5 x 5 x 5 mesh is 
used for the 8-atom zinc-blende InN and GaN unit cells 
(finer meshes resulted in total energy changes of less than 
0.1 meV/atom). Table ITTT1 shows that the relaxed lat- 
tice parameters for GaN and InN in both the wurtzite 
and zinc-blende structures are in good agreement with 
other theoretical and experimental work. The wurtzite 
GaN (InN) structures have lower total energies than zinc- 
blende GaN (InN) by 13 meV/cation (18 meV/cation). 
This is in good agreement with the calculations of Groose 
and Neugebauer, where the wurtzite GaN (InN) is lower 
in total energy than the zinc-blende GaN (InN) by 
11 meV/cation (19 meV/cation). 

We determine the equilibrium structure by minimizing 
the total energy of the system with respect to the lattice 
parameters and the positions of each of the 32 ions within 
the SQS cell. The formation enthalpy of the alloy is 

AH Y = S r (In x Gai_ x N)-a; J E; wz (InN)-(l-a;) J Bwz(GaN), 

(1) 

where x — 0, 1/4, 1/2, 3/4, or 1, and where Y represents 
either wurtzite (WZ) or zinc-blende (ZB). The wurtzite 
crystals are used as references for both the wurtzite and 
zinc-blende cases, such that the two can be easily com- 
pared. The formation enthalpy results are shown in 
Fig. □ We fit this data to the formi 6 . AH = (a + 
/3x)x(l — x) (i.e., beyond regular solution theory). For 
the wurtzite structure, we obtain a = 0.4050 cV/cation 
and /3 = —0.1117 cV/cation. For zinc-blende, we find 
a = 0.4223 cV/cation and (3 = -0.1088 eV/cation. 

Consistent with the results of Saito and Arakawai 6 - (see 
Fig. ^) , we find that all of the curves are skewed slightly 
to the left. Focusing on the wurtzite results, we see that 
the valence force field (VFF) predictions for the heat 
of formation 1 ^ are significantly smaller than those ob- 
tained here (by as much as 15.4 meV/cation at x = 1/2). 
That the heat of formation obtained from the the VFF 
is smaller than the first principles predictions was first 
noted by Grosse and Neugebauer— (by 14.2 meV/cation 
in their case). 



V. PHONON CALCULATIONS 

The phonon contribution to the free energy was de- 
termined within the harmonic approximation using the 
supercell force-constant method^L 4 ^* 4 ^^ to calculate the 
phonon density of states. In this approach, each atom 
is displaced and the forces acting on each of the other 
(static) atoms are calculated and used to determine the 
matrix of force constants. We then calculate the dynam- 
ical matrix for several reciprocal lattice vectors q and 
diagonalized these to find the phonon eigen-frequencies 
uj and vectors. Integrating the phonon eigen-frequencies 
over the Brillouin zone yields the phonon density of states 
gioS). At atmospheric pressure, the difference between 
the Gibbs and Helmholtz free energies is negligible for 
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the condensed phase. Therefore, we calculate the phonon 
free energy G v for an iV-atom supereell as a function of 
temperature T as 



G V {T) = Nk B T / duj g{uj) log 



2sinh 



Tiuj 
2k B T 



(2) 

where g(uf) is normalized according to J Q duj g(uj) = 3. 

The displacements of the atoms used to calculate the 
force constants must be small enough such that the 
system behaves harmonically yet large enough that the 
forces are sufficiently large as to be reliable. In our imple- 
mentation of the vibrational thermodynamics code, we 
employ a displacement of 0.03 A— Because of the low 
symmetry inherent to the random solid solution, each 
atom in the SQS unit cell is displaced individually. For 
the wurtzite structure, we use a supereell of 64 atoms 
(twice as large as the 16 cation SQS) to reduce the effect 
of periodic images for the determination of force con- 
stants. For GaN and InN in the wurtzite structure, we 
use a supereell consists of 3 x 3 x 3 four-atom wurtzite 
primitive cell. For the zinc-blende structure, we use 64- 
atom supercclls for all x. The LO/TO splitting, associ- 
ated with long range dipole-dipole interactions affects a 
small region near the center of the Brillouin zone, was 
not considered here (for computational efficiency). This 
should introduce only very small errors in the free energy 
of formation of the allojiS. It is possible to account for 
the LO/TO splitting using density functional perturba- 
tion theory (DFPT)^ 3 - 5 !. 

The calculated phonon densities of states of both 
wurtzite and zinc-blende In^Gai-^N are shown in Fig. [21 
for several alloy compositions. The density of states shifts 
to lower frequencies as the indium concentration is in- 
creased, as expected since (1) In atoms have a larger mass 
than Ga atoms and (2) InN has a lower— bulk modulus 
(1.39 Mbar) than GaN (2.02 Mbar). These two effects 
are of the same magnitude and together explain most of 
the shift in the vibrational frequencies in going from GaN 
to InN. Alloying InN with Ga (or vice versa) decreases 
the magnitudes of the peaks in the DOS while broadens 
the frequency range. The broadening of the phonon spec- 
tra on alloying is associated with disorder reducing the 
coherence of the phonon modes. 

We calculate the vibrational free energy difference 
AG Vi y (Y is ZB or WZ) as 

AG V , Y = G v . Y (ln x Ga 1 - x N)-xG V:W z{lnN)-(l-x)G ViWZ 

(3) 

where the same wurtzite reference frame is used. Fig- 
ure |3 shows that AG„,y is a nearly linear function of 
temperature for all compositions. For (nearly) all com- 
positions, AG„.y(i) for the wurtzite structure lies be- 
low those for the zinc-blende. The vibrational contribu- 
tion to the entropy of formation (AS V = —dAG v /dT) 
is larger for the alloys (0 < x < 1) than for x = 
or 1 (the negative AS^.zb at some compositions is an 
artifact of using wurtzite as the reference structure). 
We note that AG V for wurtzite In^Gas^N is lower 



than that for In^Gai^N. However, the AG V for zinc- 
blende In^Gaa^N is higher than that for zinc-blende 
In^Ga^N. We can fit AG v y{x) to the same type of 
interpolation function as used for the enthalpy of forma- 
tion, where the parameters a and (5 in the fit are func- 
tions of temperature. 

According to Eq. [3] AG V (x) is the difference in the 
vibrational contribution to the free energy between the 
solid solution alloy and the terminal compounds. This 
difference can be traced to the difference in phonon 
density of states (Fig. [2 and Eq. [2J between the alloy 
5in x Ga (1 _ x) N(w) and the composition weighted average 
of the density of states of GaN and InN; i.e., g av (co) = 



xgi nN (uj) + (1 - x)g GaN (uj). Comparing gi^ 



and <?av(w) in Fig. O we see that <? a v(w) is on balance 
shifted to larger lu at each x (except x = or 1, of course). 
The density of states at large frequencies u> makes a larger 
(more positive) contribution to the vibrational free en- 
ergy than do those at small frequencies. This explains 
why lattice vibrations lead to negative values of AG V . 
We return to this point below. 



VI. PSEUDO-BINARY PHASE DIAGRAM 

The free energy of formation AG of the In^Gai-^N is 
given by 



AG = AH - TAS + AG V . 



(4) 



As described above, AH and AG„ in this equation are 
known from fitting the DFT and phonon calculations 
to particular functional forms. AS corresponds to the 
entropy of mixing and, in the present analysis, is de- 
scribed in the Bragg- Williams approximation as AS — 
— ksixlogx + (1 — x) log(l — x)) (per cation). As is well- 
known, the Bragg- Williams approximation can be severe. 
However, the resulting errors are important, predomi- 
nantly in frustrated systems and systems where the phase 
diagram has a complex topology^. Neither of these sit- 
uations apply to the In^Gai-^N system. The Bragg- 
Williams approximation was also applied in many of the 
earlier calculations of the In x Gai_ x N phase diagram dis- 
cussed abovo 15 ! 16 ! 17 ! 18 . 

We employ the same wurtzite reference state for the 
unmixed GaN and InN in order to make a meaningful 
comparison of the free energies of In x Gai_a;N. Figure 0] 
(Gaffi^ws AG for both the wurtzite and zinc-blende struc- 
tures. The wurtzite structure has a lower free energy 
than the zinc-blende structure at T = 1000 K by ap- 
proximately 0.15 eV/cation. Further numerical investi- 
gation reveals that wurtzite structure has a smaller free 
energy than the zinc-blende structure at all temperatures 
and composition x investigated. This suggests that the 
wurtzite structure is always more stable than the zinc- 
blende structure in this alloy system, in agreement with 
experimental observations. 

The equilibrium solubility limits (the binodal curve) 
as a function of temperature arc calculated using a com- 
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mon tangent approach^, as shown in Fig. 0] for 1000 K. 
Figure [3] shows that the critical temperature below which 
phase separation will occur at some composition, T c , is 
2132 K for the wurtzite structure in the absence of vi- 
brational contributions to the free energy. Zinc-blende 
In-cGai-^N shows T c = 2231 K when vibrational effects 
are excluded. These values of T c (for both the wurtzite 
and zinc-blende structures) are higher than other VFF 
rcsults^^^ that neglect vibrations. These differences 
are attributable to the fact that the formation enthalpy 
determined from our DFT calculations are larger than 
those obtained within the VFF (e.g., see Fig. [TJ. Our 
predicted T c = 2231 K (for zinc-blende structure without 
vibrations) is considerably higher than those obtained by 
Teles et alM (T c = 1295 K). This is probably attributable 
to the relatively small (8-atom) supercells that they em- 
ployed. We note that another DFT calculation^ 3 ., with 
larger clusters (64 atoms/cell), predicts a slightly higher 
T c = 1400 K. The discrepancy that still remains is likely 
the result of the better description of the random solid 
solution alloy used in our analysis. 

Inclusion of the lattice vibrational effect reduces (see 
Fig. El T c for the wurtzite In x Gai_ x N from 2132 K to 
1654 K, a decrease of 29%. For zinc-blende structure, 
T c is reduced from 2231 K to 1771 K. The change in 
the critical temperature upon inclusion of the vibrational 
contributions suggests that omission of this term may be 
one of the most significant errors commonly made in the 
calculation of phase diagrams of compound semiconduc- 
tor alloys. Since this contribution is largely entropic and 
is larger in the disordered phase, it will generally lead 
to a decrease in the stability of the ordered phase rel- 
ative to the disordered phase and, therefore, to a lower 
T c . The origin of this effect can be traced to the broad- 
ening of the vibrational density of state with increasing 
disorder, as shown in Fig.|2 Inclusion of the vibrational 
contribution to the free energy also makes the wurtzite 
In x Gai_ x N phase diagram more symmetric with a peak 
at x = 0.47 with vibrations, as compared to x = 0.39 
without. The maximum in the binodal in the zinc-blende 
diagram also shifts closer to x = 0.5 with inclusion of vi- 
brational effects (although in this case the effect is weaker 
than in the wurtzite case) . We also note that inclusion of 
vibrations does not significantly increase the solubilities 
of In in GaN or Ga in InN, except above approximately 
1200 K. 

Another effect of vibrations is to decrease the width 
of the 2-phase field in the phase diagram at all temper- 
atures (relative to the phase diagram that does not in- 
clude the vibration effect). A common goal in the growth 
of In^Gai-zN alloys for microelectronic applications is 
to increase the range of compositions that are achiev- 
able without phase separation (i.e., increase the solubil- 
ity limits). Examination of the phase diagram in Fig. [3] 
shows that this could be accomplished simply by increas- 
ing the growth temperature. This is not always possible. 
An alternative approach is to increase the entropy of the 
disordered alloy relative to that of the phase separated 



material. This could be accomplished through alloying 
(quaternary additions). Addition of a new component 
increases the configurational entropy of the alloy. The 
present results suggests that alloying could also be used 
to increase the vibrational contribution to the stability of 
the disordered phase. This would be most effective with 
the addition of heavier elements, which shift the vibra- 
tional frequencies to lower values and therefore increase 
the width of the vibrational density of state (see Fig.[5J). 
Of course, the new alloy component should not increase 
the enthalpy of formation of the disordered phase. 

Figure shows a comparison of the wurtzite and zinc- 
blende In x Gai_ x N phase diagrams (including the effects 
of lattice vibrations). The wurtzite structure has a lower 
T c and is more symmetric (with respect to composition) 
than is the zinc-blende structure. These differences can 
be important in cases where the zinc-blende structure is 
stabilized through heteroepitaxyi2*i^ or other means. 

A direct comparison of the predicted phase diagram 
with an experimental phase diagram cannot be made 
since the latter has never been determined. However, an 
experimental study^ of Ino.09Gao.91N showed that this 
alloy does not phase separate at the growth tempera- 
ture of 923 K. According to the phase diagram in Fig. 
this composition and temperature lies between the bin- 
odal and spinodal. This suggests that the solid solution 
is metastable, with respect to phase separation. Addi- 
tional experimental measurements^ show that phase sep- 
aration occurs in samples grown at both Ino.37Gao.63N 
and Ino.35Gao.65N and 998 K. This is also consistent with 
the predicted phase diagram — these growth conditions 
are within the spinodal region, where the disordered al- 
loy is unstable. Since these materials were produced 
by MOCVD, there is some uncertainty as to whether 
the phases observed correspond to the bulk equilibrium 
phases at this temperature. In such cases, the phases 
that are observed may be influenced/inherited from the 
surface structure and by hetero-cpitaxial strains. How- 
ever, a solid solution cannot be stable within the spinodal 
region of the phase diagram, independent of how the ma- 
terial was produced (slow diffusional kinetics may limit 
the rate at which such an instability grows). 

The phase diagrams presented here constructed based 
upon several assumptions. These include a harmonic de- 
scription of the phonons in a relatively small unit cell, a 
Bragg- Williams approximation for the entropy of mixing 
and that the solid solution has no short-range order. All 
approximations employed in this work can be improved 
by using larger unit cells, a quasi-harmonic model for the 
vibrational contributions to the free energy, more accu- 
rate treatment of the entropy of mixing (e.g., using a 
cluster variational approach^) and a self-consistent de- 
termination of the short-range order. Despite the nature 
of the approximations made, the predicted results repre- 
sent the most rigorous determination of the In x Gai_ x N 
phase diagram to date. 
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VII. SUMMARY 

We have employed accurate density-functional calcula- 
tions to study the thermodynamic properties of both the 
zinc-blende and wurtzite structures of In^Gax^N. The 
disordered alloy was modeled using the quasi-random 
structure (SQS) approach to ensure that the order pa- 
rameter accurately describes a random solid solution (the 
detailed SQS structure information is provided for the 
16-cation unit cells for both the zinc-blende and wurtzite 
structures). The vibrational contribution to the free en- 
ergy has been calculated using first-principles calcula- 
tions of the phonon spectra within the harmonic approx- 
imation. Inclusion of these effects leads to a 29% (26%) 
reduction in the critical temperature for phase separation 
in the wurtzite (zinc-blende) structure. We found that 
the wurtzite structure is thermodynamically more stable 
than zinc-blende structure at all temperatures and alloy 
compositions. The predicted miscibility gap and spin- 
odal curve are consistent with experimental observations 
on MOCVD In x Gai_ x N. Nonetheless, care must be exer- 
cised in interpreting the phases observed following such 
growth processes in terms of an equilibrium phase dia- 
gram. 

The present results suggest that inclusion of the vibra- 
tional contribution to the free energy of formation of com- 
pound semiconductor alloys can have a profound effect on 



the critical temperature for the order-disorder transition 
on the cation sublattice. This issue is particularly im- 
portant for the increasingly complex compound semicon- 
ductors of interest today. Not only does the vibrational 
entropy shift the critical temperature, it also modifies 
the solubility limits of each component and does so in a 
non-symmetric manner. The present results also point to 
the importance of determining all of the contributions to 
the free energy self-consistently; that is using the same 
energetics description for relaxing the atomic structures, 
calculating enthalpies of formation and vibrational con- 
tributions to the free energy. While the present results 
represents a significant advance in compound semicon- 
ductor alloy phase diagram calculations, further improve- 
ments should include the effects of short range order in 
the disordered phase and better approximations to the 
entropy of mixing (see e.g., Ref . l57j) . 
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TABLE I: The SQS wurtzite (hep) structure with 16 cations 
and 16 anions for A x f3i- X N alloys where A and B represent 
In and Ga cations, respectively, with x = 1/4 (or x — 3/4) 
and x — 1/2. The lattice vectors ai, a 2 , and a3 of the struc- 
tures are linear combinations of the lattice vectors (i.e., a, 
b, and c) of a hexagonal primitive cell. In Cartesian coordi- 
nates a = (ao, 0, 0), b = (— ao/2, aovS/2, 0) and c = (0, 0, c ). 
The positions of the A cations are represented in the table 
as A(a,/3, 7) corresponding to an A cation at the position 
(aai + /3a 2 + 7a 3 )/48. 



x = 1/4 


x = 1/2 


ai = 2a - 2b - 2c 
a 2 = a + 2b 
a 3 = a + b + c 

A(39,44,6), A(18,40,36) 
A(24, 16,0), A(27,44,30) 
B(0, 16,0), B(3,44, 30) 
B(6,40, 12), B(9,20,42) 
B(42,40,36), B(45,20, 18) 
B(12, 16,24), B(15,44, 6) 
B(36, 16,24), B(21,20, 18) 
B(30,40, 12), B(33,20,42) 


ai = 2a + 2b - 2c 
a 2 — 2a + b — c 
a 3 = 2a - 2b - 2c 
A(12, 16,28), A(30,8,26) 
A(24, 16, 16), A(42,8, 14) 
A(42,8,38), A(30,8,2) 
A(36, 16,28), A(0, 16, 16) 
B(0, 16,40), B(18,8,38) 
B(36, 16,4), B(6,8,2) 
B(24, 16,40), B(12, 16,4) 
B(6,8,26), B(18,8, 14) 



TABLE III: Lattice parameters of the wurtzite (WZ) and 
zinc-blende (ZB) structures for GaN and InN from the present 
study and other studies, u describes the displacement of the 
N (0001) planes from the Ga (0001) planes in units of c as 
described, in detail, in e.g., Ref. l58l 



a(A) c(A) u 



GaN (WZ) Present calc. 


3.145 


5.121 


0.377 




CalcS 


3.196 


5.206 


0.376 




Calc^ 


3.20 


5.22 


0.376 




CalciSi 


3.133 


5.108 


0.377 




CalcS 


3.162 


5.142 


0.377 




ExptiSi 


3.190 


5.189 


0.377 


InN (WZ) 


Present calc. 


3.518 


5.690 


0.379 




Calc* 


3.545 


5.761 


0.376 




Calc* 


3.48 


5.64 


0.378 




Calc* 


3.501 


5.669 


0.3784 




Expt£ 


3.544 


5.718 




GaN (ZB) 


Present calc. 


4.443 








Calc* 


4.447 






InN (ZB) 


Present calc. 


4.964 








Calc* 


4.932 







TABLE II: The 16 cation SQS zinc-blende (fee) structure for 
A^Bi-^iV alloys. The lattice vectors ai, a 2 , and a 3 are re- 
ported in terms of the lattice vectors of the primitive rhom- 
bohedral cell: a = a (0, 1/2, 1/2), b = a (l/2, 0, 1/2) and 
c = ao(l/2, 1/2, 0). The positions of the A cations are repre- 
sented in the table as A(a,f3, 7) corresponding to an A cation 
at the position (aai + /3a 2 + 7a 3 )/16. 



x = 1/4 

ai = 2a + b - 2c 
a2 = a + b + 2c 
a 3 = 2a - b + 2c 
A(ll,12,7), A(4,0,4) 
A(15,12, 11), A(9,4, 13) 
B(0,0,0), B(3, 12, 15) 
B(13,4,l), B(6,8, 14) 
B(10,8,2), B(14, 8,6) 
B(l,4,5), B(8,0,8) 
B(2,8,10), B(5,4,9) 
B(12,0,12), B(7,12,3) 



x = 1/2 

ai = 2a — 2c 
a 2 = 2a + 2b 
a 3 = 2a - 2b 

A(8,4,4), A(0,8,8) 
A(8,0,8), A(0,4,12) 
A(8,0,0), A(0,4,4) 
A(8,12,4), A(0,0,8) 
B(0,0,0), B(8, 12,12) 
B(0,12,4), B(8,8,0) 
B(0, 12, 12), B(8,8,8) 
B(0,8,0), B(8,4, 12) 
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FIG. 1: (Color online) The formation enthalpy of wurtzite 
(WZ) and zinc-blende (ZB) In^Gai^^N as a function of com- 
position x. The solid curve is a fit to the calculated values 
(open circles) with (a + f3x)x(l — x). The valence force field 
results of Saito and Arakawa— are also included for compar- 
ison. 
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FIG. 2: (Color online) Phonon density of states (DOS) of (a) 
wurtzite and (b) zinc-blende In^Gai-^N correspond to x = 0, 
1/4, 1/2, 3/4, and 1, which have been normalized according to 
g{v)du — 3, where v = lv/2-k. An arbitrary shift is added 
to the curves for clarity. The black and grey curves refer to 
the solid solution system gin ;E Ga (1 _ x) N(^') and an average of 
the InN and GaN phases g a v(w), respectively. 
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FIG. 3: (Color online) Vibrational free energy difference AG V 
as a function of temperature, for the wurtzite (solid lines) and 
zinc-blende (dashed lines) structures. 
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FIG. 4: (Color online) The common tangent intersects the 
wurtzite AG curve at xi — 0.024 and X2 = 0.95, which gives 
the binodal curve in Fig.|S| Lattice vibrational effect has been 
taken into account. A representative case of T = 1000 K is 
used. 




FIG. 5: (Color online) The pseudo-binary diagram (in 
temperature-cation concentration space) for (a) wurtzite and 
(b) zinc-blende In^Gai-zN with and without the inclusion 
of lattice vibration effects. The solid lines correspond to the 
equilibrium binodal curves and the dashed lines represent the 
spinodal curves. 
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FIG. 6: (Color online) A comparison of the pseudo-binary 
diagrams of zinc-blende and wurtzite In^Gai-^N (including 
the effects of lattice vibrations). The solid lines correspond to 
the equilibrium binodal curves and the dashed lines represent 
the spinodal curves. The points represented by the filled and 
open circles represent the experimental data^ for which phase 
separation was and was not observed, respectively. 



